#
library(parallel)
library(tictoc)
library(HeadR)
library(SQUAREM)
library(fixest)
library(xtable)
library(latex2exp)
library(Rcpp)
setwd("~/Dropbox/BLP_REStat/Replication")
rm(list=ls())
source("Rfiles/CF_BLPdata_functions.R")
source("Rfiles/CF_MMX_functions.R")
source("Rfiles/CF_MMX_MLOG_functions.R")
#
Rcpp::sourceCpp("Cppfiles/crossDelta.cpp")
#parallel processing set up for doParallel, doRNG
MAXIT <- 500
seedn <- 777 # 666
sequences <- "normal" # alternatives: normal, halton or "sobol"
#
n.hh <- 100000 # careful, this is slow, but works
U0 <- 1  # needed in prob.mat(), but not prob.mat.namedxvars()
n.reps <- 1 # for richsubs (no need to do a lot if a lot of consumers.)
#
###
# Setting up the stage ====
###
#
setwd("~/Dropbox/BLP_REStat/Replication")
source("Rfiles/CF_BLPdata_getdata.R") # get original data
#
BLP[,rank := frank(-s),by=year]
#View(BLP[year==1990,.(m,p,s,rank)][order(-s)])
year.gph <- 1990
year.monte <- year.gph
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix stuff
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
# initialization phase
set.seed(seedn)
data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=3) #gen exog data with v=3 for initialization
BLP.output <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective="richsubs") 
BLP.output.m <- BLP.output[m==j]
BLP.meanownelas <- BLP.output.m[, mean(own_elas_m)] 
BLP.priceadjust <- BLP.output.m[, mean(p_m * (1-s_m))] 
#
###
# Calibrating alpha  ====
###
#
ysd.values <- seq(0,2,by=0.02) # the values for ysd.
alpha_interval <- c(0.5,100)
#
# Function
alpha.search <- function(ysd){
  #  
  alpha.loop <- function(alpha.guess,ysd){
    vsim <- 3
    y.sd <<- ysd
    alpha.data <- alpha
    alpha <<- alpha.guess
    set.seed(seedn)
    data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=vsim) 
    BLP.output <- xi.BLP.mlog(dat=data.exog.mm,verbose=FALSE) 
    eps.E.mean <- superelas_mean(dat=data.exog.mm,xi.dat=BLP.output)
    alpha <<- alpha.data
    return(eps.E.mean$eps.mean)
  }
  #find the right alpha to match the average elasticity
  ufn <- function(x) alpha.loop(ysd=ysd,alpha.guess=x)-BLP.meanownelas
  alpha.root <- uniroot(ufn,interval=alpha_interval)$root
  return(data.table(ysd,alpha.root))
}
#alpha.search(1.72)
#
tic()
DT.alpha <- mclapply(ysd.values,alpha.search,mc.cores = detectCores()-1)  
DT.alpha <- rbindlist(DT.alpha)
toc()
#resume
setorder(DT.alpha,ysd)
DT.alpha
#
###
# Mean superelas approach ====
###
#
# Function
ysd.loopE <- function(ysd.loop){
  vsim <- 3
  y.sd <<- ysd.loop
  alpha.data <- alpha
  alpha <<- DT.alpha[ysd==ysd.loop]$alpha.root
  set.seed(seedn)
  data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=vsim) 
  BLP.output <- xi.BLP.mlog(dat=data.exog.mm,verbose=FALSE) 
  eps.E.mean <- superelas_mean(dat=data.exog.mm,xi.dat=BLP.output)
  alpha <<- alpha.data
  return(data.table(ysd.loop,eps.E.mean))
}
#
# Run the loop
tic()
DT <- mclapply(ysd.values,ysd.loopE,mc.cores = detectCores()-1)  
DT <- rbindlist(DT)
toc()
setorder(DT,ysd.loop)
DT
saveRDS(DT,"Data/RDS/superE_pte.rds")
DT = readRDS("Data/RDS/superE_pte.rds")
#
brilliantblue =rgb(0,117,197,maxColorValue = 255)
#
#margins to use everywhere
pltmgn = c(4,4,1,1)
pdf(file="Graphs/CF_BLPdata/superE_pte_ysd101.pdf",height=4,width=4.5)
par(mar=pltmgn)#bottom, left, top and right
DT[,plot(ysd.loop,pte.mean,col="black",type="l",lwd=2,ylim=c(-0.55,1.3),xaxs="i",
         xlab="s.d. of log income",ylab=TeX("mean elasticity"))]
DT[,lines(ysd.loop,E.mean,col="grey50",type="l",lwd=2)]
abline(h=1,lty="dashed")
abline(h=0,lty="dashed")
abline(v=1.72,lty="dotted")
abline(v=0.9,lty="dotted")
mtext(side=3,at=0.9,"PTE=1, E=0")
mtext(side=3,at=1.72,"BLP sd")
text(0.62,0.25,"Super Elas.",col="grey50",pos=2)
text(1.24,1.33,"Pass-through",col="black",pos=1)
dev.off()
#
